Relaxation optimized transfer of spin order in Ising spin chains 
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In this manuscript, we present relaxation optimized methods for transfer of bilinear spin correla- 
tions along Ising spin chains. These relaxation optimized methods can be used as a building block 
for transfer of polarization between distant spins on a spin chain. Compared to standard techniques, 
significant reduction in relaxation losses is achieved by these optimized methods when transverse 
relaxation rates are much larger than the longitudinal relaxation rates and comparable to couplings 
between spins. We derive an upper bound on the efficiency of transfer of spin order along a chain 
of spins in the presence of relaxation and show that this bound can be approached by relaxation 
optimized pulse sequences presented in the paper. 
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I. INTRODUCTION 

Relaxation (dissipation and decoherence) is a charac- 
teristic feature of open quantum systems. In practice, 
relaxation results in loss of signal and information and ul- 
timately limits the range of applications. Recent work in 
optimal control of spin dynamics in the presence of relax- 
ation has shown that these losses can be significantly re- 
duced by exploiting the structure of relaxation 0,000]. 
This has resulted in significant improvement in sensitiv- 
ity of many well established experiments in high resolu- 
tion nuclear magnetic resonance (NMR) spectroscopy. In 
particular, by use of optimal control methods, analytical 
bounds have been achieved on the maximum polarization 
or coherence that can be transferred between coupled 
spins in the presence of very general decoherence mecha- 
nisms. In this paper, we look at the more general problem 
of transfer of coherence or polarization between distant 
spins on an Ising spin chain in the presence of relax- 
ation. This problem is ubiquitous in multi-dimensional 
NMR spectroscopy, where polarization is transferred be- 
tween distant spins on a chain of coupled spins. Spin 
(or pseudo spin) chains also appear in many proposed 
quantum information processing architectures |3, |6| . 

The system that we study in this paper is a linear chain 
of n weakly interacting spins 1/2 placed in a static exter- 
nal magnetic field in the z direction (NMR experimental 
setup), with Ising type couplings of equal strength be- 
tween nearest neighbors, see Fig. ^ The free evolution 
Hamiltonian of the system has the form 



27rJ hzl{i+i)z, 

i=l 



where u>i is the Larmor frequency of spin i and J is the 
strength of the coupling between the spins. In a suitably 
chosen (multiple) rotating frame, which rotates with each 
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FIG. 1: The system that we study in this paper is a linear 
chain of n weakly interacting spins 1/2 with Ising coupling 
between next neighbors. The coupling constant is the same 
for all pairs of connected spins. 



spin at its resonance (Larmor) frequency, the free evolu- 
tion Hamiltonian simplifies to 



n-i 
i=i 



(1) 



Motivated by NMR spectroscopy of large molecules 
in solution, we assume that the relaxation rates of the 
longitudinal operators with components only in the z 
direction, like Ii z , 2Ii Z Iu+i-\ z is negligible compared 
with relaxation rates for transverse operators like Ii X , 
2IizI(i+i)y 0- The transverse relaxation is modeled by 
the Lindbladian with the general form 



L(P) 
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In liquid state NMR spectroscopy, the two terms of 
the Lindblad operator model the relaxation mechanism 
caused by chemical shift anisotropy and dipole-dipole in- 
teraction respectively 0. Here we neglect any interfer- 
ence effects between these two relaxation mechanisms 0] . 
Relaxation rates , bij depend on various physical pa- 
rameters, such as the gyromagnetic ratios of the spins, 
the internuclear distance, the correlation time of molecu- 
lar tumbling etc. We define the net transverse relaxation 
rate for spin i as k l a — fli+J^j Without loss of gen- 
erality in the subsequent analysis, we assume that k l a are 
equal and we denote this common transverse relaxation 
rate by k. 



The time evolution (in the rotating frame) of the spin 
system density matrix p is given by the master equation 



P = -i[H,p] + L(p) 



(2) 



where H = H c + H r f and H r f is the control Hamiltonian. 
In the NMR context, the available controls are the com- 
ponents of the transverse radio-frequency (RF) magnetic 
field. It is assumed that the resonance frequencies of the 
spins are well separated, so that each spin can be selec- 
tively excited (addressed) by an appropriate choice of the 
components of the RF field at its resonance frequency. 

Consider now the problem of optimizing the polariza- 
tion transfer 



Iu 



(3) 



along the linear spin chain shown in Fig. ^ in the 
presence of the relaxation mechanisms mentioned above. 
This problem can be stated as follows: Find the opti- 
mal transverse RF magnetic field in the control Hamil- 
tonian H r f such that starting from p(0) = I\ z and 
evolving under Eq. J5J) the target expectation value 
(Inz) — kr{p(t)I nz } is maximized. 
To fix ideas, we analyze the case when n = 3 



Iu 



hz 



(4) 



This transfer is achieved conventionally using INEPT like 
pulse sequences 0,^,^3- Under the conventional trans- 
fer method, the initial state of the system I\ z evolves 
through the following stages 



hz — * 2-Zi z /2 



2I 2z hz 



'3z- 



(5) 



In the first stage of the transfer, I\ z — > 2Ii z I 2zi spin 3 
is decoupled from the chain using standard decoupling 
methods [12| and the initial polarization I\ z on spin 1 is 
rotated by RF field (an appropriate 7r/2 pulse) to coher- 
ence l\ x , which then evolves under coupling Hamiltonian 
2I\ z Iiz to 21\ y I 2z . When the expectation value (2I\ y Ii z ) 
is maximized, another 7r/2 pulse is used to rotate 2I\ v l 2z 
to 2/i 2 /22- This is the INEPT pulse sequence. The next 
stage, 2I lz I 2z — * 2I 2z I 3z , is the so-called spin order trans- 
fer. Fig. |3 shows the population inversion corresponding 
to this transfer. By a suitable ir/2 rotation of spin 2, 
the density operator 2I lz I 2z is transformed to 2I lz I 2x , 
which then evolves to 2I 2x I 3z . When the expectation 
value (2I 2x I 3z ) is maximized, another ir/2 pulse is used 
to rotate 2I 2x I 3z to the spin order 2I 2z I 3z . This is the 
Concatenated INEPT (CINEPT) pulse sequence. The fi- 
nal stage transfer, 2I 2z I 3z — > I 3z , is similar to that in the 
first stage and is accomplished by the INEPT pulse se- 
quence. The efficiency of these transfers is limited by the 
decay of transverse operators, due to the phenomenon of 
relaxation. 

In our recent work on relaxation optimized control of 
coupled spin dynamics , we showed that the efficiency 
of the first and the last step (the two INEPT stages) 
in Eq. (jSJ) can be significantly improved by controlling 
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FIG. 2: Energy level diagram for the spin order transfer 
2Ii z l2z — » 2l2 Z l3z- The dark circles represent population ex- 
cess. State a corresponds to spin up and state /3 to spin down. 



precisely the way in which magnetization is transferred 
from longitudinal operators to transverse operators. In 
other words, instead of using 7r/2 (hard) pulses to rotate 
longitudinal to transverse operators (and the inverse), we 
can exploit the fact that longitudinal operators are long 
lived by making these rotations gradually, saving this way 
magnetization. This transfer strategy is called relaxation 
optimized pulse element (ROPE). 

In this article we derive relaxation optimized pulse 
sequences for the intermediate transfer (the CINEPT 
stage), 



2/i z /2 



2hJ, 



(6) 



This relaxation optimized transfer of spin order can then 
be used as a building block for the polarization transfer 
© through the scheme 



Iiz ~ * 2I\ Z I 2 



2I,J: 



2z*3z 



21, 



In 



(7) 



II. THE OPTIMAL CONTROL PROBLEM AND 
AN UPPER BOUND FOR THE EFFICIENCY 

In this section, we formulate the problem of trans- 
fer in Eq. © as a problem of optimal control and 
derive an upper bound on the transfer efficiency. To 
simplify notation, we introduce the following symbols 
for the expectation values of operators that play a part 
in the transfer. Let z\ = (2I\ z I 2z ), x\ = (2Ii z I 2x ), 
y 2 = (V2{2I lz I 2y I 3z + I 2y /2)), x 3 = -(2I 2x I 3z ) and 
23 = (2I 2z I 3z ). As a control variable we use the trans- 
verse RF magnetic field, pointing say in the y direction 
(in the rotating frame), so H r f = ui y (t)I 2y . Note that 
u! y (t) is the component of the field in the rotating frame, 
so it is actually the envelope of the RF field. The carrier 
frequency of the RF field is the resonance frequency of 
spin 2. Using Eq. @, we find that the evolution of the 
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system, in time units of 1/ , is given by 
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(8) 



where fi„(i) = w v (t)/(irJV2) and £ = k/{jV2). The 
initial condition is (zi, Xi,y2, x 3 , 23) — (1, 0, 0, 0, 0). 

The efficiency of the conventional method (CINEPT) 
for transfer Z\ — > Z3, can be easily found. At t — 0, z\ is 
transferred to X\ by application of a (tt/2)^ pulse on spin 
2. Couplings evolve xi to j/2 which further evolves to X3. 
As a function of time, 2:3 (t) = e - ^' sin 2 (t/v / 2)- This is 
maximized for t m = v / 2cot~ 1 (^/v / 2)- At i = i m a second 
{tt/2) v pulse is applied on spin 2, rotating the maximum 
value from 23 to 23. This value is the efficiency 77c"/ of 
the conventional method 



Vci 



= exp(-^V2cot- 1 (^/v / 2))sin 2 (cot- 1 (,e/v / 2)) . (9) 



A better efficiency can be achieved if we store magne- 
tization in the decoherence free longitudinal operators Zi 
while the system is evolving. This is done by rotating z% 
to xi gradually, instead of using tt/2 hard pulses. This 
is the physical concept behind the relaxation optimized 
transfer strategy. For the specific transfer examined in 
this article, we first find an upper bound for the maxi- 
mum efficiency and in the next section we calculate nu- 
merically the magnetic field £l y (t) that approaches this 
bound. 

In order to derive the upper bound, we use an aug- 
mented system instead of the original one iJHJ. The aug- 
mentation is done in two steps. First, we suppose that 
we can rotate z\ to X\ and X3 to Z3 independently using 
two different controls, say Qi(t) and SI3 (i) , instead of the 
common control Q y (t). Next, we provide 7/2 with a relax- 
ation free partner zi and with a control Q2 (i) which can 
rotate yi to z%. The augmented system is 



" Zl ' 




" 








-Oi 





- 




" Zl ' 


Zl 
















-fia 







Z2 


Z3 























Z3 


±1 




ill 








-e 


-1 







Xl 


2/2 







a 2 





t 


-£ 


-f 




2/2 


. £3 . 




. 





<h 





1 


-€ - 




. X3 . 



(10) 

Observe that system (fTU)) reduces to system © for fii = 
— O3 = fly and CI2 — 0. Thus, if we know the maximum 
achievable value of z 3 starting from (1,0,0,0,0,0) and 
evolving under system l|10|l , then this is an upper bound 
for the maximum achievable value of Z3 with evolution 
described by the original system JSJ. 

and 



Let n 
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+ zf, 
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\J x\ + z\. Using fli(t) we can control the angles 9i, 
shown in Fig. [21 independently. If we assume that the 
control can be done arbitrarily fast as compared to the 
evolution of couplings or relaxation rates then we can 
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FIG. 3: Auxiliary variables rt. 



think of 9i as control variables. The equations for the 
evolution of r, are 
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(11) 



where the new control parameters are Ui(t) = cos(9i(t)). 
The goal is to find the largest achievable value of T3 start- 
ing from (ri,r2,r3) = (1,0,0) by appropriate choice of 
Ui{t). This problem can be solved analytically. The opti- 
mal solution is characterized by maintaining vanishingly 
small values of dr2/dt, i.e., (uiri — u$rs)/u2T2 = £ and 
by U3T3/uiri = k, where 



(12) 



The maximum achievable value of T3 is also k. We prove 
it in the following. 

Using variables p t = rf, m, = uin/ \Jy^ (u.n) 2 , and 

dr/dt = (uiTi) 2 , equation (jll|) can be re-written as 



d_ 

dr 



where 



Pi( T ) 
M T ) 



A : 



diag{Arn{r)m r {t)) 



-i -1 
1 -i -1 

1 -£ 



(13) 



(14) 



m T = (mi,TO2,m3) and diag(X) represents the vector 
containing diagonal entries of the square matrix X. The 
goal is to find the controls TOj(r) m f( T ) — 1) and the 
largest achievable value of P3 starting from {pi,P2,P3) — 
(1,0,0). Eq. O implies that 



diag(A / m(r)m T \r)dr) 
Jo 



Let M — m(T)m T (r)dr . Note that M is a symmet- 
ric, positive semidefinite matrix. By definition Pi(r) > 
for all t. At final time T, we must have Pi(T) = and 
P2{T) = 0, as any nonzero value of Pi(T) or piiT) can 



'MT) 




Mo) 






Mo) 






Mo) 



4 



be partly transferred to and the final value of Ps{T) 
further increased. Since (pi(0),p 2 (0),P3(0)) = (1,0,0), 
it implies that M should be such that [AM)n = — 1, 
(AM)%2 — and (AM)^ is maximized over all positive 
scmidefinite M satisfying the above constraints. This 
problem is a special case of a semidcfmitc programming 
problem |13| . If the symmetric part of matrix A is neg- 
ative definite, as in our case l|14|) . then it can be shown 
that the optimal solution M to the above stated semi- 
definite programming problem is a rank one matrix |l4j , 
i.e., M = mm T for some constant m and therefore the 
ratio Uiri/ujTj in (|llfl is constant throughout. The con- 
dition (AM)22 = 0, implies (u\ri — u^r^)/^ = U2T2- Sub- 
stituting for U2T2 in JUJ, we obtain 
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(15) 



We now simply need to maximize the gain in r% to loss 
in r%, i.e., the ratio r^/i—ri). This yields u^r^/uiri = k, 
with k given in O- The corresponding maximum effi- 
ciency for transfer rj — > is also ft. This is the maximum 
efficiency for transfer z\ — > Z3 under the augmented sys- 
tem l|10|). and thus an upper bound for the efficiency of 
the same transfer under the original system 



TABLE I: For various values of £ £ [0, 1], the optimal values of 
A, a and the corresponding efficiency are shown. We present 
also for comparison the efficiency achieved by the steepest 
descent method. 



I 


A 


a 


Gaussian Pulse 


Steepest Descent 


1.00 


1.11 


1.30 


0.2510 


0.2512 


0.95 


1.09 


1.32 


0.2661 


0.2662 


0.90 


1.07 


1.34 


0.2824 


0.2825 


0.85 


1.05 


1.36 


0.3000 


0.3001 


0.80 


1.03 


1.38 


0.3190 


0.3191 


0.75 


1.02 


1.39 


0.3396 


0.3397 


0.70 


1.00 


1.41 


0.3619 


0.3620 


0.65 


0.98 


1.43 


0.3861 


0.3863 


0.60 


0.97 


1.44 


0.4124 


0.4126 


0.55 


0.96 


1.44 


0.4410 


0.4413 


0.50 


0.95 


1.44 


0.4721 


0.4726 


0.45 


0.94 


1.45 


0.5060 


0.5067 


0.40 


0.93 


1.46 


0.5428 


0.5439 


0.35 


0.92 


1.46 


0.5830 


0.5846 


0.30 


0.91 


1.46 


0.6270 


0.6292 


0.25 


0.90 


1.47 


0.6750 


0.6780 


0.20 


0.89 


1.48 


0.7277 


0.7315 


0.15 


0.88 


1.48 


0.7855 


0.7900 


0.10 


0.85 


1.52 


0.8494 


0.8536 


0.05 


0.79 


1.60 


0.9203 


0.9232 


0.00 


0.73 


1.71 


0.9999 


1.0000 



III. NUMERICAL CALCULATION OF THE 
OPTIMAL RF FIELD AND DISCUSSION 

Having established an analytical upper bound (|12fl for 
the efficiency, we now try to find numerically a RF field 
£l y (t) that approaches this bound, for each value of the 
parameter £. We emphasize that in this section the orig- 
inal system JSJ is employed. 

At first, we use a numerical optimization method based 
on a steepest descent algorithm. For the application of 
the method we use a finite time window T. For values of 
normalized relaxation £ in (0—1), a time interval T = 
10 (normalized time units) is enough. For larger values 
of £ we can use even shorter T. The optimal RF field 
Q y {t) that we find with this method, for various £, is 
shown in Fig. 0] Note that as £ increases, the optimal 
pulse becomes shorter in time and acquires a larger peak 
value. The reason for this is that for larger £ the transfer 
z\ — > Z3 should be done faster, in order to reduce the time 
spent in the transverse plane and hence the relaxation 
losses. Now observe that the optimal pulse shape can be 
very well approximated by a Gaussian profile of the form 



fl y (t) — Aexp 



t-T/2 
V2a 



(16) 



with A, a appropriately chosen. As a result, the effi- 
ciency that we find using the appropriate Gaussian pulse 
is very close to that we find using the original pulse. This 
suggests that instead of using the initial numerical opti- 
mization method, we can use Gaussian pulses of the form 



Hlfty . optimized with respect to A and a for each value of 
£. The optimal A, a are found by numerical simulations. 
For each £ we simulate the equations of system (jHJ) with 
Fly(t) given by Eq. (|16fl . for many values of A and a. 
We choose those values that give the maximum 23 (T). 
In table I, we show the optimal A, a for various values 
£ £ [0, 1]. We also show the corresponding efficiency, as 
well as the efficiency achieved by the initial numerical 
optimization method. Observe how close lie these two 
groups of values. The choice of the Gaussian shape is 
indeed successful. 

Fig. 13 shows the efficiency of the conventional method 
(CINEPT), i.e., i] C i from Eq. JjJ), the efficiency of our 
method (SPORTS ROPE, SPin ORder TranSfer with Re- 
laxation Optimized Pulse Element) , and the upper bound 
k from Eq. I|12l) , for the values of relaxation parameter £ 
shown in table I. Note that for large £ (large relaxation 
rates), SPORTS ROPE gives a significant improvement 
over CINEPT. Also note that it approaches fairly well 
the upper bound. 

Using the Gaussian pulse shape we can get a quantita- 
tive impression of the robustness of SPORTS ROPE. In 
Fig. [S] we give a gray-scale topographic plot of the effi- 
ciency (2:3 (T) with T = 10) as a function of A and a for 
£ = 1. The maximum value can be found from table I and 
is 0.2510. The white region corresponds to values > 0.24, 
while the black region to values < r\ci{£, = 1) = 0.1727. 
The intermediate gray regions correspond to values be- 
tween these two limits. Obviously, SPORTS ROPE is 
quite robust. 

In Fig. Ufa), we plot the time evolution of the various 
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FIG. 4: Optimal pulse (dashed line) calculated using a numerical optimization method based on a steepest descent algorithm, 
for various values of the normalized relaxation parameter £. The Gaussian pulse (solid line) approximates very well the optimal 
pulse shape and gives a similar efficiency. This suggests that instead of using the initial numerical optimization method, we 
can use Gaussian pulses of the form I16H . optimized with respect to A and a for each value of £. 



transfer functions (expectation values of operators) that 
participate in the transfer z\ — > z 3 , when the optimal 
Gaussian pulse for £ = 1, shown in Fig. Hfc), i s applied to 
system ©. Observe the gradual building of the interme- 
diate variables x±, y% and X3. Note that dr^jdt — 1/2 7^ 0. 
There is no contradiction with the optimality condition 
dr 2 /dt = derived in section [n] during the calculation of 



the upper bound, since this condition refers to the aug- 
mented system (I1U|I and not the original one JSJ used 
here. In Fig. Efb) we plot the angle #3 = tan -1 (Z3/X3) 
of the vector r% with the x axis, as a function of time. 
Observe that initially r 3 is parallel to x axis (8 3 = 0), but 
under the action of the Gaussian pulse is rotated gradu- 
ally to z axis (#3 = 7r/2). This gradual rotation of r 3 (as 
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FIG. 5: Efficiency for the conventional method (CINEPT), 
Eq. and for our method (SPORTS ROPE), for the values 
of £ shown in table I. The upper bound (1121 for the efficiency 
is also shown. 



well as of n) is a characteristic feature of the SPORTS 
ROPE transfer scheme. 

We remark that for the general transfer I\ z — > 
I nz , more than one intermediate steps 2Iu^x)zhz ~ * 
2IizI(i+i) z are necessary. Since the equations that de- 
scribe the i th transfer are the same as (jHJ, we just need 
to apply the same Gaussian pulse but centered, in the 
frequency domain, at the resonance frequency of spin i. 
In this sequence of Gaussian pulses we should add at the 
beginning and at the end the optimal pulses for the first 
and the final step, respectively, see Eq. (J7J). These pulses 
can be calculated using the theory presented in Fi- 
nally, note that the same scheme can be used for the 
coherence transfer I± a — > I n p, where a,/3 can be x or y. 
We just need to add the initial and final tt/2 pulses that 



accomplish the rotations I\ a — > I\ z , I n 



IV. CONCLUSION 
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In this paper, we derived an upper bound on the ef- 
ficiency of spin order transfer along an Ising spin chain, 
in the presence of relaxation, and calculated numerically 
relaxation optimized pulse sequences approaching this 
bound. Using these methods, a significant reduction in 
relaxation losses is achieved, compared to standard tech- 
niques, when transverse relaxation rates are much larger 
than the longitudinal relaxation rates and comparable 
to couplings between spins. These relaxation optimized 
methods can be used as a building block for transfer of 
polarization or coherence between distant spins on a spin 
chain. This problem is ubiquitous in multi-dimensional 
NMR spectroscopy and is also interesting in the context 
of quantum information processing. 



FIG. 6: Gray-scale topographic plot of the efficiency 
as a function of the amplitude A and the standard 
deviation a of the Gaussian pulse for f = 1. The max- 
imum value can be found from table I and is 0.2510. 
The white region corresponds to values > 0.24, while 
the black region to values < r)ci(£=l) = 0.1727. The 
intermediate gray regions correspond to the intervals 
[0.23, 0.24), [0.22, 0.23), [0.21, 0.22), [0.20, 0.21), [0.1727, 0.20) 
(from white to black). 
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of the vector 7-3 with the X cLXlS, ELS £L function of time. Observe that initially r-j, is parallel to x axis (63 = 0), but under the 
action of the Gaussian pulse is rotated gradually to z axis (63 — tt/2). 
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